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We present a computationally tractable scheme of time-dependent transport phenomena within 
open-boundary time-dependent density-functional-theory. Within this approach all the response 
properties of a system are determined from the time-propagation of the set of "occupied" Kohn- 
Sham orbitals under the influence of the external bias. This central idea is combined with an open- 
boundary description of the geometry of the system that is divided into three regions: left/right 
leads and the device region ( "real simulation region" ) . We have derived a general scheme to extract 
the set of initial states in the device region that will be propagated in time with proper trans- 
parent boundary-condition at the device/lead interface. This is possible due to a new modified 
Crank-Nicholson algorithm that allows an efficient time-propagation of open quantum systems. We 
illustrate the method in one-dimensional model systems as a first step towards a full first-principles 
implementation. In particular we show how a stationary current develops in the system indepen- 
dent of the transient-current history upon application of the bias. The present work is ideally 
suited to study ac transport and photon-induced charge-injection. Although the implementation 
has been done assuming clamped ions, we discuss how it can be extended to include dissipation due 
to electron-phonon coupling through the combined simulation of the electron-ion dynamics as well 
as electron-electron correlations. 



PACS numbers: 72.10.-d, 73.23.-b, 73.63.-b 
I. INTRODUCTION 

During the last decades, the size of electronic cir- 
cuits has continuously been reduced. Today, systems 
like quantum wires and quantum dots are routinely pro- 
duced on the nanometer scale. Recently, the seemingly 
ultimate limit of minituarization has been achieved by 
several experimental groups who were able to place sin- 
gle molecules between two macroscopic electrodesi^ To 
describe transport properties on such a small scale, a 
quantum theory of transport is required^* 

A cornerstone of such a theory is the Landauer- 
Biittiker formalism 5,6 which provides a method to com- 
pute the steady-state current of non-interacting electrons 
for meso- or nanoscopic systems connecting two (or more) 
macroscopic electrodes. 

Alternatively, the technique of nonequilibrium Green 
functions (NEG)i£ has been used to tackle quantum 
transport. Studies using the NEG approach typically 
use tight-binding- like model Hamiltonians to describe the 
combined system electrodes plus "device" . A well known 
scheme is the one introduced by Caroli et o/Aif. In the 
remote past the left and right electrodes are disconnected 
and in equilibrium at two different chemical potentials; 
the conducting part of the Hamiltonian is switched on 
adiabatically and eventually a steady-state develops. 
Within this contacting approach also time-dependent 
transport phenomena have been investigated^ Caroli et 
al. discussed non-interacting systems only. Their ap- 



proach has later been extended to account for short-range 
electron-electron interaction and for interaction with vi- 
brations in the device region^ An excellent overview of 
the field can be found in the book by Haug and Jauhoi^ 
and in Ref. 0. Despite its appeal, the above scheme has 
limitations since the time-dependent perturbation is the 
tunneling Hamiltonian, connecting the electrodes to the 
device, rather than the external electric field. 

Cini proposed another scheme^ also based on NEG. 
Here, the system electrodes plus "device" is connected 
and in equilibrium in the remote past. The time- 
dependent perturbation is the external scalar potential. 
It has been showni 5 " that for non-interacting systems the 
contacting approach and the Cini approach yield the 
same current in the long-time limit and that in the dc 
case the steady state current does not depend on the his- 
tory of the applied potential. Moreover, the Cini scheme 
is well suited for a density functional extension since the 
electrons are driven out of equilibrium by a local poten- 
tial rather than by a non-local one (see below). 

With recent experimental progress to place single 
molecules as devices between macroscopic electrodes 
there also has been considerable activity to describe 
transport through these systems on an ab initio level. 
Most approaches are based on a self-consistency proce- 
dure first proposed by Lang*i& In this steady-state ap- 
proach based on density functional theory (DFT), ex- 
change and correlation is approximated by the static 
Kohn-Sham (KS) potential and the charge density is 
obtained self-consistently in the presence of the steady 
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current. However, the original justification involved 
subtle points such as different Fermi levels deep inside 
the left and right electrodes and the implicit reference 
of non-local perturbations such as tunneling Hamiltoni- 
ans within a DFT framework. (For a detailed discus- 
sion we refer to Ref. lia.) The steady-state DFT ap- 
proach has been further developecUZiiSiiLSli and the re- 
sults have been most useful for understanding the quali- 
tative behavior of measured current-voltage characteris- 
tics. Quantitatively, however, the theoretical I-V curves 
often differ from the experimental ones by several orders 
of magnitude^ Several explanations are possible for such 
a mismatch: models are not sufficiently refined, parasitic 
effects in measurements have been underestimated, the 
characteristics of the molecule-contact interface are not 
well understood and difficult to address given their atom- 
istic complexity. Adding to the theoretical reason for 
this discrepancy is the fact that the transmission func- 
tions computed from static DFT have resonances at the 
non-interacting Kohn-Sham excitation energies which in 
general do not coincide with the true excitation energies. 
Furthermore, different exchange-correlation functionals 
lead to DFT-currents that vary by more than an order 
of magnitude^ 

Excitation energies of interacting systems are acces- 
sible via time-dependent (TD) DFT. 23,24 In this theory, 
the time-dependent density of an interacting system mov- 
ing in an external, time-dependent local potential can be 
calculated via a fictitious system of non-interacting elec- 
trons moving in a local, effective time-dependent poten- 
tial. Therefore this theory is in principle well suited for 
the treatment of noncquilibrium transport problems^ A 
basic issue is that most exchange-correlation functionals 
have been derived under equilibrium conditions and their 
application to non-equilibrium problems should be ana- 
lyzed in more detail. However, this is beyond the scope 
of the present work. 

Before a TDDFT calculation of transport can be tack- 
led, a number of technical problems have to be addressed. 
In particular, one needs a practical scheme for the prop- 
agation of the time-dependent Schrodinger equation for 
an infinitely large system. Of course, since one can in 
practice only deal with finite systems this can only be 
achieved by applying the correct boundary conditions. 
The problem of so-called "transparent boundary condi- 
tions" for the time-dependent Schrodinger equation has 
been attacked by many authors. For a recent overview, 
the reader is referred to Ref. I2H 

In this paper we present a propagation scheme which 
is particularly designed to be used for the calculation of 
time-dependent transport problems. In Section [H] wc 
combine the Cini scheme with TDDFT and we develop 
a general formalism based on the propagation of Kohn- 
Sham orbitals in open systems. In Section IIIII we will 
address the question of how to obtain the correct initial 
states for the propagation. An algorithm for the time- 
evolution of open systems is proposed in Section llVI It 
is based on a modified version of the Crank-Nicholson 
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FIG. 1: Schematic diagram of the system described in the 
main text for the calculation of charge transport through the 
central constriction. 



algorithm. Section [V] describes some details of our nu- 
merical implementation and Section IVII gives results for 
several one-dimensional model systems. We draw our 
conclusions in Section IVni 



II. GENERAL FORMULATION 

We consider an electrode-junction-electrode system 
which is initially in equilibrium (t < 0). The system is 
contacted and no current flows through the junction, the 
charge density of the electrodes being perfectly balanced, 
see Fig. ^ Therefore, the system initially is in its ground 
state which, due to the Hohenberg-Kohn theorem^ is 
a functional of the density. This density can then be 
computed in the usual way by n(r, 0) = J2 occ |-0 S (r, 0)| 2 
where the sum is over the occupied Kohn-Sham or- 
bitals , 0s(r, 0), i.e., the eigenfunctions of the Kohn-Sham 
Hamiltonian H(0) with eigenergy below the Fermi en- 
ergy. Here and in the following we use boldface notation 
to denote operators in one-electron Hilbert space. 

To observe a current we drive the system out of equi- 
librium by exposing the electrons to an external elec- 
tric potential (bias). Without loss of generality we 
will assume that this external bias vanishes for times 
t < 0. According to the Runge-Gross theorem^ 3 ., the 
time-dependent density can be calculated by evolving the 
KS orbitals according to the KS equation of TDDFT 
iips(t) = H(t)ip s (t) where H(t) is the time-dependent 
KS Hamiltonian. Thus, n(r, t) — ^ occ | V's ( r ) t)\ 2 an d the 
continuity equation is J^n(r,t) = —V • jKs( r :i), where 
jxsM) = -E OM MVS(r>t)W»(r,t)] is the KS cur- 
rent density. Due to the Runge-Gross theorem and the 
continuity equation one can deduce that the longitudinal 
part of the KS current density equals the longitudinal 
part of the true current density. This need not be true 
for the transverse part of the current density. However, 
the transverse part of the current density does not con- 
tribute to the total current which can then be calculated 
by a surface integral 



Is(jt) 



N 



d<rA-Im[V>;(r,t)W,(r,t)], (1) 



where n is the unit vector perpendicular to the surface 
element da and the surface S is perpendicular to the lon- 
gitudinal geometry of our system. In order to propagate 
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the KS orbitals we need to solve the Schrodinger equation 
for a macroscopic and non-periodic system. This goal 
is hopeless unless we know the dynamics of the remote 
parts of the system. We restrict ourselves to metallic 
electrodes. Then, the external potential and the distur- 
bance introduced by the device region are screened deep 
inside the electrodes. As the system size increases, the 
remote parts are less disturbed by the junction and the 
density inside the electrodes approaches the equilibrium 
bulk-density. Thus, the macroscopic size of the electrodes 
leads to an enormous simplification since the initial-state 
self-consistency is not disturbed far away from the con- 
striction. 

It is convenient to partition the system into three main 
regions: a central region C consisting of the junction 
and a few atomic layers of the left and right electrodes 
and two regions L, R which describe the left and right 
bulk electrodes. Only the central device region C will 
be treated explicitly. Our scheme accounts for the full 
dynamical screening in the central region. It can be fur- 
ther refined by taking into account screening effects also 
deeper in the electrodes at the level of linear response, 
with a limited increase in numerical efforts. (These ef- 
fects might be of importance in the initial transient phase 
where long-range plasma oscillations in the electrodes 
may occur.) According to the above partitioning, the 
original KS Hamiltonian can be written as a 3 x 3 block 
matrix, and the Schrodinger equation reads 

. d 

where ip a is the projected wave-function onto the region 
a = L,R,C. We can solve the diffential equation for tpL 
and ipR by introducing the retarded Green function g a for 
electrode a, which satisfies {z^ — H aa {t)\ g a (t, t') = 
5{t — t'), a — L, R, with boundary conditions g a (t + , t) 
i and g a (t, t + ) = 0. Then, we have for a = L,R 

7jj a (t) =z<7 Q (t,O)0 Q (O)+ / dt'g a (t,t')H aC ip c (t'). (3) 

Jo 

Using Eq. ©, the equation for tpc can be written as 

i§-^c(t) = H CC {t)i>c{t) + J dt'E(i,t')W) 

+i J2 H Ca g a (t,0)M0), (4) 

a=L,R 

where £ = J2 a =L,R H c a g a H aC is the self-energy which 
accounts for the hopping in and out of region C. Thus, 
for any given KS orbital we can evolve its projection 
onto the central region by solving Eq. in region 
C. Eq. (@J has also been derived elsewhere (for static 
Hamiltonians)^ To summarize, all the complexity of 
the infinite electrode-junction-electrode system shown in 
Fig. n has been reduced to the solution of an open 
quantum-mechanical system (the central region C) with 
proper time-dependent boundary conditions. 
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III. COMPUTATION OF EXTENDED 
EIGENSTATES 

Eq. Q of the previous Section is the central equation 
of our approach to time-dependent transport. It is a re- 
formulation of the original time-dependent Schrodinger 
equation J5J of the full system in terms of an equation 
for the central (device) region only. The coupling to 
the leads is taken into account by the lead Green func- 
tions g a , a = L,R. Eq. has the structure of a time- 
dependent Schrodinger equation with two extra terms: 
the term involving the self energy £ which we will also 
denote as the memory integral since it involves the wave- 
function in the central region at previous times during 
the propagation. 29 The second term describes the injec- 
tion of particles induced by a non- vanishing projection of 
the initial wave- function onto the leads. 

Eq. 10} is first order in time, therefore we need to 
specify an initial state which is to be propagated. We 
want to study the time evolution of systems perturbed 
out of their equilibrium ground state. Of course, the 
ground state of our noninteracting system is the Slater 
determinant of the occupied eigenstates of the full, ex- 
tended Hamiltonian in equilibrium, H(t < 0) = H s . 
The practical question then arises how one can obtain 
these eigenstates without having to deal explicitly with 
the extended Hamiltonian. Note that, unlike in a bulk 
solid, the translational symmetry is broken in the present 
device situation. 

In the present Section we propose a solution of this 
problem which is based on the partitioning approach 
used in many steady-state transport calculations (see, 
e.g., Ref. l3l|) . The retarded Green function of the static 
Hamiltonian in the energy domain is determined by 



[(E + i V )l-H s ]g(E) = l 



(5) 



The Green function Q(E) of the full system can be writ- 
ten in the same block structure as the Hamiltonian 



0(E) = 



Gll(E) Qlc(E) Qlr(E) 
Qcl{E) g C c(E) Gcr(E) 
Qrl{E) Q rc{E) Qrr(E) 



(6) 



Eq. (O can be solved for the block of the Green function 
referring only to the central region 



Qcc(E) 



1 



(E + i v )l c - H s cc - Z a =L,R HS c a 9 a (E)H s aC 

with the retarded Green function of lead a 

1 



g a (E) = 



{E + 11])^ - H & 



(7) 



(8) 



and the unit matrix 1 Q in region a. This Green func- 
tion enters as a central ingredient into the Fisher-Lee 
relation^ for the calculation of the transmission func- 
tion. Through the coupling to the leads it provides for 
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level broadening of the isolated central part, but it also 
contains information on the eigenstates of the extended 
system. 

In order to illustrate the central idea of our method to 
extract the extended eigenstates from the Green function 
we consider H s to be the discretized form of a continuous 
Hamiltonian H s (r). The continous Green function and 
the discretized one for uniform lattice spacing Ax are 
connected by 



g{vi,v 3 ,E) = g(E) l3 /(Aa 



(9) 



where N = 1,2,3 is the number of spatial dimensions 
of the problem. We choose the convention that a single- 
particle orbital ipEj of the Hamiltonian H s is uniquely 
specified by its eigenenergy E and a label j for the cIe 
degenerate orbitals of this energy. Using the Lehmann 
representation and assuming that H s is invariant under 
time-reversal, the imaginary part of Q is 

-~Im [g{v,v',E)}=Y,6{E-E')^ Elj {v)r Elj {r') . 

E> 3 = 1 

(10) 

Multiplying Eq. ((1111) by ipEi(*'), integrating over r', using 
the orthogonality of the single particle states and dividing 
by the density of states D(E) = J2e> $(E ~~ E')d E > we 
obtain 



1 



nD(E) 
with 



d N r Im [g(r,r',E)]ijEi(r')= 7 (E)ijEi(r) 

(11) 

(12) 



E' 



Eq. (|llfl has the structure of an eigenvalue equation 
where the energy eigenstate ipEi is also an eigenstate 
of the integral operator — Im [0(r, r', E)]/ (irD(E)) with 
eigenvalue ^(E). For a given energy E, this integral op- 
erator has cIe different degenerate eigenstates. 

We note that Eq. (|10f) is valid for all points in space, 
in particular also for both r and r' representing points 
in the central region. In this case we know the Green 
function Qcc through Eqs. J7J) and 10. Below we show 
that the eigenfunctions of Im[^cc] can be expressed as 
linear combination of the ipEi- Let us consider the matrix 
formed by the elements 



b m i = / d N r ip% m (r)ipEi(T) 
Jc 



(13) 



where the integration is over the central region only. This 
matrix is Hermitian and can be diagonalized, i.e., 



^ brnldjl — Xj 



(14) 



with Xj real. Next we compute the matrix elements of 
the Green function with respect to the functions 



aEj(r) =^2aji{E)ipEi(r) ■ 



(15) 



i=i 



After a straightforward manipulation one finds 

-- [d N r fdV a H (r)Im[g(r,r',£)]a Bi (r') 
w Jc Jc 

= 5 j i^ y £5(E-E') (16) 

E' 

which shows explicitly that the functions a£j(r) diago- 
nalizc Im [Qcc] m the central region and that the eigen- 
values are positive. Since any linear combination of de- 
generate eigenstates is again an eigenstate, diagonalizing 
Im [Qcc{E)\ gives us one set of linearly independent, de- 
generate eigenstates of energy E. 

In our practical implementation described in more de- 
tail in Section we diagonalize 



TrD c (E) 



Im [g CC {E)} 



where 



D C {E) = — Tr {Im [Qcc(E)}} 

IT 



(17) 



(18) 



is the total density of states in the central region. If we 
use N g grid points to describe the central region, the di- 
agonalization in principle gives N g eigenvectors but only 
a few have the physical meaning of extended eigenstates 
at this energy. It is, however, very easy to identify the 
physical states by looking at the eigenvalues: only few 
eigenvalues (for the simple examples we studied either 
one or two) are nonvanishing and they always add up to 
unity. The corresponding states are the physical ones. 
All the other eigenvalues are zero (or numerically close 
to zero) and the corresponding states have no physical 
meaning. 

The procedure described above gives the correct ex- 
tended eigenstates only up to a normalization factor. 
When diagonalizing Eq. 117|) with typical library rou- 
tines one obtains eigenvectors which are normalized to 
the central region. Physically this might be incorrect. 
Therefore, the normalization has to be fixed separately. 
In the example of Section^we fixed the norm by match- 
ing the wavefunction for the central region to the known 
form (and normalization) of the wavefunction in the 
macroscopic leads. 

It should be emphasized that the procedure described 
here for the extraction of eigenstates of the extended sys- 
tem from Qcc{E) only works in practice if E is in the 
continuous part of the spectrum due to the sharp peak of 
the delta function in the discrete part of the spectrum. 

Eigenstates in the discrete part of the spectrum can be 
found by considering the original Schrodinger equation 
for the full system: 



H s i/j = Eijj. 



(19) 
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Using again the block structure of the Hamiltonian this 
can be transformed into an effective Schrodinger equa- 
tion for an energy- dependent Hamiltonian for the central 
region only: 



H cc - ^2 H c a g a {E)H aC ] ipc = Eip C - (20) 

a=L,R 

This equation has solutions only for certain values of E 
which are the discrete eigenenergies of the full Hamilto- 
nian H s . Therefore, one can find these states by iter- 
ation. We have succesfully tested this idea for systems 
where the analytic solutions are known. However, since 
the main focus of the present work is transport where the 
continuum states are the essential physical ingredient, we 
will not deal with the states in the discrete spectrum for 
the remainder of this paper. Those states might play a 
role in the description of charge-accumulation in molecu- 
lar transport, as, e.g., in Coulomb blockade phenomena. 



IV. ALGORITHM FOR TIME EVOLUTION 

In order to calculate the longitudinal current in an 
electrode-junction-electrode system we need to propa- 
gate the Kohn-Sham orbitals. The main difficulty stems 
from the macroscopic size of the electrodes whose remote 
parts, ultimately, are taken infinitely far away from the 
central, explicitly treated, scattering region C . 

The problem can be solved by imposing transparent 
boundary conditions^ at the electrode-junction inter- 
faces. Efficient algorithms have been proposed for wave- 
packets initially localized in the scattering region and 



for Hamiltonians constant in time. In this Section we 
propose an algorithm well suited for delocalized initial 
states, as well as for localized ones, evolving with a time- 
dependent Hamiltonian. 

Let H{t) be the time-dependent KS Hamiltonian. We 
partition H (t) as in Section^ The explicitly treated re- 
gion C includes the first few atomic layers of the left and 
right electrodes. The boundaries of this region are chosen 
in such a way that the density outside C is accurately de- 
scribed by an equilibrium bulk density. It is convenient to 
write H aa (t), with a = L, R, as the sum of a term H s aa 
which is constant in time and another term U a (t) which 
is explicitly time-dependent, H aa (t) — H s aa + U a (t). In 
configuration space U a (t) is diagonal at any time t since 
the KS potential is local in space. Furthermore, the diag- 
onal elements U a (r,t) are spatially constant for metallic 
electrodes. Thus, U a (t) = U a (t)l a and U L (t) - U R (t) 
is the total potential drop across the central region. We 
write H(t) = H(t) + U(t) with 



and 



H(t) 



U(t) = 



H) 



Hi 







LL ±± LC 

Hcl Hcc(t) Hen, 
H RC H s rr 



u L (t)u o 


U R (t)l R 



(21) 



(22) 



In this way, the only term in H(t) that depends on t 
is Hcc(t)- For any given initial state ip(0) = r/>'°) we 
calculate tp(t m — to At) = ij;( m > by using a generalized 
form of the Cayley method 



(l + i&ET ; ) 



- (mKl+ifC/ (m) 



(m) 



(m+1) 



(1-^)1-4^^ 



(m) 



(23) 



withi* (m) = \[H% 



i)+H(t r 



i[t/(Wi) + 



U(t m )] and S = At/2. It should be noted that our 
propagator is norm conserving (unitary) and accurate to 
second-order in 8, as is the Cayley propagator.? 3 Denot- 
ing by ip a the projected wave function onto the region 
a = R,L,C, we find from Eq. 



Here, H 



central region: 
i8H s LL )-^H 



\-i8H 



(m) 
off 



1 + iSH 



(m) 
off 



+ 5 (to) - M (m) . 



(24) 



(to) 
eff 



is the effective Hamiltonian of the 



LC 



H 



(m) 
CC 



;[Hcc(tm+l 



i8H CR {l + i8H 

Hcc(tm)} 



i8H CL {\ + 
1 Hrc: with 



RR) ±J -RC 

The source term 



g(m) describes the injection of density into the region 
C, while the memory term is responsible for the 

hopping in and out of the region C. In terms of the prop- 
agator for the uncontacted and undisturbed a electrode 



9a 



1 - iSHl 
1 + iSm 



the source term can be written as 



S (r, 



2iS 



(m,0) 



1 + iSH, 



(m) 

eff a=L,R 



Am) 



-Hca 



(25) 



l9 a \ 



(0) 



i + ism 



(26) 
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with onto the left and right electrode tpa^ vanishes and — 

. . for any m, as it should be. The memory term is more 

«W = 1 ~ % f° m) and A^ fe ) = n^)] 2 . (27) complicated and reads 

For a wave packet initially localized in C the projection 



M< m ) = -- 



1 + iSH 



(m) 

eff a=L,R k=0 



(m) O) 



[Qa 



(m—k) 



Qi" 



Wc 



(28) 



where 



Q(m) rr [ffaj ii 

1 + l6H nr , 



(29) 



The quantities Q„ depend on the geometry of the sys- 
tem and are independent of the initial state 

Below we propose a recursive scheme to calculate the 
Q^^'s for those system geometries having semiperiodic 
electrodes along the longitudinal direction, see Fig. |2J In 
this case H s aa has a tridiagonal block form 



and are given by 



a {m) _ y 

Ha — v a 



l9 a \ 



1 + iSH n 



(32) 



i,i 



where the subscript (1,1) denotes the first diagonal block 
of the matrix in the square brackets. We introduce the 
generating matrix function 



1 



x + iy5H no 



' a i 



(33) 



h a V a 
V a h a V a 

V a h a 



(30) 



where h a describes a convenient cell and V a is the hop- 
ping Hamiltonian between two nearest neighbor cells. 
Without loss of generality we assume that both h a and 
V a are square matrices of dimension N a x N a . Taking 
into account that the central region contains the first few 
cells of the left and right electrodes, the matrix has 
the following structure 



(m) 





















(31) 



The q^'s are square matrices of dimension N a x N a 
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FIG. 2: Schematic sketch of an electrode-junction-electrode 
system with semiperiodic electrodes. 



which can also be expressed in terms of continued matrix 
fractions (see Appendix^) 



i 



x + iySha + y 2 5 2 q a {x,y) 



V, 



The q^'s can be obtained from 



_(m) = _±_ 



d_ d_ 

dx dy 



(34) 



(35) 



x=y=l 



From Eqs. (|35|l and (|34|) one can build up a recur- 
sive scheme. Let us define p~ 1 (x,y) = x + iydh a + 

y 2 S 2 q a (x,y) and p^ = + ^baO> J/)U=s/=i- 

Then, by definition, 



Ha — y aP a y c 



(36) 



Using the identity — 



ml L dx dy 



] m Pa(x,y)p a 1 (x,y) = o, 



one finds 

(1 + l5h a )Pa m) = (1 - i5h a )Pa m ' 1) 

m 

-5 2 J>« + 2 ^- 1} + q { t 2) )P [ a- k) (37) 



fc=0 



with p^ — qa^ = for m < 0. Once q^ has been 
obtained by solving Eq. (|34|) with x = y = 1, we can 



calculate p a 



(<>) 



iSh a 

J 1 



S 2 qa^} 1 . Afterwards, we 



can use Eq. I|37|l with q a — V a Pa V a to calculate p a 
and hence qi 1 and so on and so forth. 



This concludes the description of our algorithm for the 
propagation of the time-dependent Schrodinger equation 
for extended systems. It is worth mentioning an addi- 
tional complication here which arises for the propagation 
of a time-dependent Kohn-Sham equation. This com- 
plication stems from the fact that in order to compute 
ipc +1 ^ at time step m + 1 one needs to know the time- 
dependent KS potential at the same time step which, via 
the Hartree and exchange-correlation potentials, depends 

on the yet unknown orbitals %1)^q 1+1 \ Of course, the solu- 
tion is to use a predictor-corrector approach: in the first 

step one approximates Hqq by H cc (t m ) , computes new 
orbitals and from those obtains an improved ap- 

proximation for Hqq . 



V. IMPLEMENTATION DETAILS 

All the methodological discussion above is general and 
can be applied to general device configurations as long as 
they can be mapped into a longitudinal-like geometry as 
described in Fig. [3 In order to demonstrate the feasibil- 
ity of the scheme described in the previous Sections we 
have implemented it for one-dimensional model systems. 
The extension to real molecular-device configurations is 
presently under development i24 We have used a simple 
three-point discretization for the second derivative 



1 



(Ax) 



(ip(x i+ i) - 2ip(xi) + ip(xi-i)) 

(38) 

with equidistant grid points Xi, i = 1, . . . , N g and spac- 
ing Ax. Within this approximation matrices of the form 
HcaMHac which are N g x N g matrices and appear, 
e.g., in Eq. Q or i|29[l . have only one nonvanishing ma- 
trix element. For a = L this is the (1,1) element, for 
a = R it is the [N g , N g ) element. 

In order to proceed we have to specify the nature of 
the leads and therefore the lead Green function. Here we 
choose the simplest case of semi-infinite, uniform leads at 
constant potential U a Q. In this case, the Green function 
(jSJ) in the energy domain can be given in closed form: 



[g a (E)h -- 

iAx 



iAx | 
— ==exp<^ i\]2E a \x k - xi\ 

V2E a I 



\J2E a 



exp <^ iy 2E a (\x k - x a0 \ + \xi - x a0 \ 



•39) 



with E a = E — U a Q. The abscissa x a o is the position of 
the interface between the lead and the device region and 
Xk — x a o ± kAx, where the plus sign applies for a = R 
and the minus sign for a = L. 

The results of the procedure for calculating extended 
eigenstates as described in Section 11111 is illustrated in 
Fig. |21 for a square potential barrier with zero potential 
in both leads. In the upper panel we have the square 
modulus of eigenstates at an energy below the barrier 





FIG. 3: Continuum states of square potential barrier at dif- 
ferent energies with leads at zero potential. Upper panel: 
eigenstates for e = 0.45 a.u., just below the barrier height of 
0.5 a.u.. Lower panel: eigenstates at e = 0.6 a.u.. 



height while in the lower panel eigenstates with energy 
higher than the barrier are shown. The states result from 
diagonalization of Eq. I|17|) . In order to obtain the nor- 
malization constant we compute the logarithmic deriva- 
tive at the boundary of the central region numerically 
and match it to the analytic form in the lead to obtain 
the phase shift 5 a : 



i__E_ 

2dx 2 



H\m\ 2 ) 



qcot(S a ) 



(40) 



where q = \f 2E a . Knowing the phase shift we can rescale 
the wavefunction such that it matches with the analytic 
form sm.(q(x — x a o) +6 a ) at the interface. Of course, this 
form of the extended states only applies for E a > but 
as long as E is in the continuous part of the spectrum, it 
is correct at least for one of the leads. This is sufficient 
to determine the normalization. The states obtained nu- 
merically with this procedure coincide with the known 
analytical results. 

We then implemented the propagation scheme pre- 
sented in the previous Section. Within our three-point 
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U=0.05 

U-0.05 steady state 
U=0.15 
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U=0.25 
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FIG. 4: Time evolution of a Gaussian wavepacket with ini- 
tial width 1.0 a.u. and initial momentum 0.5 a.u. for various 
propagation times. The transparent boundary conditions al- 
low the wavepacket to pass the propagation region without 
spurious reflections at the boundaries. 



approximation, h ai V a and q a arc lxl matrices. The 

equation for [see Eqs. (I3"4l) and becomes a sim- 
ple quadratic equation which can be solved explicitly 



(0) 



-(l + iSh a ) + y/(i + i8KY + 4((5vg 2 

25 2 



(41) 



Although the quadratic equation has two solutions, the 

above choice for qa^ is dictated by the fact that the Tay- 
lor expansions for small i5 of Eqs. (|41|l and (|34|l have to 
coincide. Using this result we then solved the iterative 
scheme to obtain the q£"^ for m > 1. 

As a first check on the propagation method we prop- 
agated a Gaussian wavepacket which, at initial time 
t = 0, is completely localized in the central device region. 
(The source terms then vanish identically). As can 
be seen in Fig. 21 the wavepacket correctly propagates 
through the boundaries without any spurious reflections. 

For the propagation of the extended initial states (the 
eigenstates of the unperturbed system) we also need to 
implement the source terms S( m \ In the following we 
assume that the left and right leads are at the same po- 
tential initially so that the analytic form of the initial 
states is in both leads given by sin(q(x — x a o) + S a ) — 
cos(5 a ) sin(q(x — x a o)) + sin(5 Q ) cos(q(x — x a o)). The 
propagation of the term proportional to sm(q(x — x a o)) 
is trivial since this is an eigenstate of the lead Hamilto- 
nian with energy e q = q 2 /2. Therefore, if - in discretized 

form - i/j r1 = (sin(qAx), sin(2qAx), . . .) T we obtain: 



H 



cn 



[9rI 



1 + iSH 



RR 



-^ (0) 
YRA 



Vr 



(1 - i6e q ) m 
(l + iSeq) m+1 



en 



(42) 



H S RR is the static 



and similarly for the left lead. Here 
part of the right-lead Hamiltonian, g R the corresponding 
Green function according to definition (|25|l and e R = 
(0, . . . , 0, 1) T is a unit vector. 



FIG. 5: Time evolution of the current for a system where 
initially the potential is zero in the leads and the propagation 
region. At t = 0, a constant bias with opposite sign in the 
left and right leads is switched on, U = Ul — —Ur (values in 
atomic units). The propagation region extends from x — —6 
to x — +6 a.u.. The Fermi energy of the initial state is 
ef = 0.3 a.u.. The current in the center of the propagation 
region is shown. 



The propagation of the part proportional to cos(q(x — 
%ao)) is more complicated since this is not an eigenstate 
of the lead Hamiltonian. We define the function R(x, y) 
from 



R(x,y)e R = H C r 



x + iy6H s RR 



^ (0) 

VR.2 



(43) 



where ip R2 — (cos(gAx), cos(2gAa;), . . .) T . Introducing 
the tridiagonal matrix 



O 



/ 1 •••\ 
1 1 ••• 
10 1- 



V 



(44) 



and using 



VrO^ = (e q - h R )^ 2 - V R (1, 0, 0, . . .f (45) 
one arrives at 



,(0) 



R(x,y) = 



1 



[Vr cos(gAx) + iy5q R (x, y)} , (46) 



x + iySe q 

where qn(x,y) is given by Eq. (|33|) . Defining 



x=y=l 



(47) 



finds 



+«E „ (1 ~f;C t (^'+4 t -" 

^ (1 + i5e„) m+1 - k V R H 

k=Q y q ' 



(48) 
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and finally 

l + idH RR 

One may proceeds along the same lines for extracting the 
left component. 

To test our implementation we have propagated eigen- 
states of the extended system. As expected, these states 
just pick up an exponential phase factor exp(— iEt) dur- 
ing the propagation. 

VI. EXAMPLES 

We are now in a position to apply our algorithm 
to the calculation of time-dependent currents in one- 
dimensional model systems. The systems are initially 
in thermodynamic equilibrium. At time t — 0, a bias is 
switched on in the electrodes. 

As a first example we considered a system where the 
electrostatic potential vanishes identically both in the left 
and right leads as well as in the central region which is 
explicitly propagated. Initially, all single particle levels 
are occupied up to the Fermi energy ep . At t = a 
constant bias is switched on in the leads and the time- 
evolution of the system is calculated. We chose the bias 



in the right lead as the negative of the bias in the left 
lead, Ur = —Ul- The current is calculated from Eq. QJ: 

/ kp dk / d \ 

27 Im (^fcOM^feMj 

r k F dk / d d \ 

= 2 /„ s Im (*a;* + **'-3i*-) < 50 » 

where the prefactor 2 comes from spin and kp — \f2e~p 
is the Fermi wavevector. 

The numerical parameters are as follows: the Fermi 
energy is ep — 0.3 a.u., the bias is Ul = — Ur = 
0.05,0.15,0.25 a.u., the central region extends from x = 
—6 to x = +6 a.u. with equidistant grid points with 
spacing Ax = 0.03 a.u.. The /c-integral in Eq. (|50fl is 
discretized with 100 fc-points which amounts to a propa- 
gation of 200 states. The time step for the propagation 
was At = 10~ 2 a.u. 

In Fig. we have plotted the current densities at x = 
as a function of time for different values of the applied 
bias. As a first feature we notice that a steady state is 
achieved, in agreement with the results of Ref. Il5l The 
corresponding steady-state current / can be calculated 
from the Landauer formula. For the present geometry 
this leads to the steady current 



I = 8e 



max([/i,(/ K ) 



duj 

2^ 



\f(w-U L )-f(u-U R )]- 



yju - Ul^/u> ~ U R 



yzi^ui + - U R \ +U l U r 



sin(/\/2a7) 



(51) 



where / is the width of the central region. From Eq. I|51|) 
with I — 12 a.u. and Ul — ~Ur, the numerical values 
for the steady-state currents are 0.0316 a.u. (Ul = 0.05 
a.u.), 0.0883 a.u. (U L = 0.15 a.u.) and 0.0828 a.u. (U L = 

0. 25 a.u.). We see that our algorithm yields the same 
answers. Second, we notice that the onset of the current 
is delayed in relation to the switching time t — 0. This is 
easily explained by the fact that the perturbation at t = 
happens in the leads only, e.g., for |a;| > 6 a.u., while we 
plot the current at x = 0. In other words, we see the delay 
time needed for the perturbation to propagate from the 
leads to the center of our device region. We also note that 
the higher the bias the more the current overshoots its 
steady-state value for small times after switching on the 
bias. Finally it is worth mentioning that increasing the 
bias not necessarily leads to a larger steady-state current. 

In the second example we studied a double square po- 
tential barrier with electrostatic potential V(x) = 0.5 
a.u. for 5 a.u. < |x| < 6 a.u. and zero otherwise. This 
time we switch on a constant bias in the left lead only, 

1. e., Ur — 0. The Fermi energy for the initial state is 
ep = 0.3 a.u.. The central region extends from x = —6 



to x = +6 a.u. with a lattice spacing of Ax = 0.03 
a.u.. Again, we use 100 different fc-values to compute the 
current and a time step of At = 10~ 2 a.u.. 

In Fig. we plot the current at x = as a function 
of time for several values of the applied bias U = Ul- 
Again, a steady state is achieved for all values of U. As 
discussed in Fig. the transient current can exceed the 
steady current; the higher the applied voltage the larger 
is this excess current and the shorter is the time when it 
reaches its maximum. Furthermore, the oscillatory evo- 
lution towards the steady current solution depends on 
the bias. For higher bias the frequency of the transient 
oscillations increases. For small bias the electrons at the 
bottom of the band are not disturbed and the transient 
process is exponentially short. On the other hand, for a 
bias close to the Fermi energy the transient process de- 
cays as a power law, due to the band edge singularity. 
As pointed out in Ref. [Til for non-interacting electrons 
the steady-state current develops by means of a pure de- 
phasing mechanism. In our examples the transient pro- 
cess occurs in a femtosecond time-scale, which is much 
shorter than the relaxation time due to electron-phonon 
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FIG. 6: Time evolution of the current through a double 
square potential barrier in response to an applied constant 
bias (given in atomic units) in the left lead. The potential 
is given by V(x) = 0.5 a.u. for 5 < \x\ < 6 a.u. and zero 
otherwise, the propagation region extends from x = —6 to 
x = +6 a.u.. The Fermi energy of the initial state is sf ~ 0.3 
a.u.. The current in the center of the structure is shown. 



interactions. 

In Fig. we plot the time evolution of the total num- 
ber of electrons in the device region for the same values 
of Ul- We see that as a result of the bias a quite sub- 
stantial amount of charge is added to the device region. 
This result has important implications when simulating 
the transport through an interacting system as the effec- 
tive (dynamical) electronic screening is modified due not 
only to the external field but also to the accumulation of 
charge state in the molecular device. This illustrates that 
linear response might not be an appropriate tool to tackle 
the dynamical response and that we will need to resort to 
a full time-propagation approach as the one of the present 
work. Here we emphasize that all our calculations are 
done without taking into account the electron-electron 
interaction. If we had done a similar calculation with 
the interaction incorporated in a time-dependent Hartree 
or time-dependent DFT framework we would expect the 
amount of excess charge to be reduced significantly as 
compared to Fig. [7| 

In Fig.[S]we show time-dependent currents for the same 
double barrier as in Fig. for two different ways of ap- 
plying the bias in the left lead: in one case the constant 
bias Uq is switched on suddenly at t = (as in Fig. EJ), in 
the other case the constant Uq is achieved with a smooth 
switching U(t) = U sm 2 (ujt) for < t < n/(2w). As ex- 
pected and in agreement with the results of Ref. ITHl the 
same steady state is achieved after the initial transient 
time. However, the transient current clearly depends on 
how the bias is switched on. 

In the final example we address the simulation of ac- 
transport. We computed the current for a single square 
potential barrier with V(x) — 0.6 for \x\ < 6 and zero 
otherwise. Here we applied a time-dependent bias of 
the form Ul{€) = Uos'm(uit) to the left lead. The right 
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FIG. 7: Time evolution of the total number of electrons in 
the region \x\ < 6 for the double square potential barrier of 
Fig.H 

lead remains on zero bias. The numerical parameters 
are: Fermi energy ep — 0.5 a.u., device region from 
x = —6 to x = +6 a.u. with lattice spacing Ax = 0.03 
a.u.. The number of fc-points is 100 and the time step 
is At = 10~ 2 a.u.. In Fig. we plot the current at 
x = as a function of time for different values of the 
parameter J7 = 0.1,0.2,0.3 a.u. The frequency was cho- 
sen as uj = 1.0 a.u. in both cases. Again, as for the 
dc-calculation discussed above, we get a transient that 
overshoots the average current flowing through the con- 
striction; again, this excess current is larger the higher 
the applied voltage. Also, after the transient we obtain 
a current through the system with the same period as 
the applied bias. Note, however, that (especially for the 
large bias), the current is not a simple harmonic as the 
applied ac bias. 



VII. SUMMARY AND OUTLOOK 

In the present work we have presented a formally rig- 
orous approach towards the description of charge trans- 
port using an open-boundary scheme within TDDFT. 
We have implemented a specific time-propagation scheme 
that incorporates transparent boundaries at the de- 
vice/lead interface in a natural way. In order to have 
a clear definition of a device region in Fig. ^ we assumed 
that an applied bias can be described by adding a spa- 
tially constant potential shift in the macroscopic part 
of the leads. This implies an effective "metallic screen- 
ing" of the constriction. The screening limits the spatial 
extent of the induced density created by the bias po- 
tential or the external field to the central region. Our 
treatment can be further refined to include dynamical 
screening deep inside the electrodes on the level of linear 
response, which might be of importance for the initial 
transient. Our time-dependent scheme allows to extract 
both ac and dc I/V device characteristics and it is ide- 
ally suited to describe external field (photon) assisted 
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FIG. 8: Time evolution of the current for a double square 
potential barrier when the bias is switched on in two different 
manners: in one case, the bias Uo is suddenly switched on at 
t — while in the other case the same bias is achieved with 
a smooth switching U(t) = Uos'm 2 (uit) for < t < 7t/(2oj). 
The parameters for the double barrier and the other numerical 
parameters are the same as the ones used in Fig. |S| Uo and 
uj given in atomic units. 



FIG. 9: Time evolution of the current for a square potential 
barrier in response to a time-dependent, harmonic bias in 
the left lead, Uh(t) = Uo sin(ttrf) for different amplitudes Uo 
(values in a.u.) and frequency lu — 1.0 a.u.. The potential is 
given by V (x) = 0.6 a.u. for \x\ < 6.0 a.u. and zero otherwise. 
The propagation region extends from x — —6 to x = +6 a.u.. 
The Fermi energy of the initial state is ef = 0.5 a.u.. The 
current at x — is shown. 



processes. 

In order to illustrate the performance of the method, 
we have implemented it for one-dimensional models and 
we have shown: i) How to extract the proper initial ex- 
tended states to be propagated, ii) How to incorporate 
perfect transparent boundaries for the time propagation, 
iii) A steady-state current is always reached upon ap- 
plication of a dc bias. The transient process occurs on 
a time-scale much shorter than the relaxation time due 
to electron-phonon interaction. In the case of systems 
without any source of dissipation it is known that the 
steady-current is independent of the history of the pro- 
cess We have explicitly demonstrated this history in- 
dependence for two different switching processes of the 
external bias. However, if we allow for dissipation ei- 
ther through electron-electron or electron-phonon inter- 
actions, the current versus voltage characteristics might 
depend on the history. For instance, hysteresis loops due 
to different transient electronic/geometrical device con- 
figurations are possible. This effect will be more dramatic 
in the case of ac-driving fields of high frequencies where 
the system may not have enough time to respond to the 
perturbation, iv) A periodic ac current is reached upon 
perturbation with a monochromatic field. 

Previous work on time-dependent quantum transport 
mainly uses the idea of Caroli>2ii£ This approach is at the 
core of the Landauer derivation and has the problem of 
using different chemical potentials in different parts of the 
system. This implies that the initial state is not a ground 
state of the entire, contacted system. Furthermore, the 
time-dependent perturbation is a tunneling Hamiltonian 
which is nonlocal in space. Therefore, it cannot be com- 
bined with TDDFT since there the time-dependent po- 
tential is local. 



Here, we have used the scheme of Ciniia which can be 
combined with TDDFT: We start the calculation from 
a well-defined thermodynamic equilibrium configuration, 
therefore the scheme is thermodynamically consistent. 
Then we apply an external potential that in general is 
time dependent. By virtue of the Runge-Gross theorem, 
the time-evolution of this quantum system is completely 
determined by the knowledge of the time-dependent den- 
sity. In the steady state regime the occupation of left and 
right moving carriers is dictated by the time-dependent 
Schrodinger equation. 

Another thermodynamically consistent scheme has 
been put forward by Kamenev and Kohn^ They used 
the microscopic quantum-kinetic formulation of conduc- 
tivity and worked with the Kubo formula— in the linear 
time-dependent Hartree regime. They consider a closed 
system (ring) where the current in the system results 
from an external driving vector potential. This approach 
also overcomes the problem of having two chemical po- 
tentials. They have shown that it is possible to recover 
the Landauer result, i.e., the universal quantum of con- 
ductance 2e 2 /h independent of the length and material. 
Since the Kamencv-Kohn approach uses a vector poten- 
tial rather than a scalar potential, time-dependent cur- 
rent DFT rather than TDDFT would be the natural 
density-functional extension. 

Most theoretical approaches to transport in molecular 
electronic devices adopt open boundary conditions and 
assume that transport is ballistic, i.e., under steady state 
conditions inelastic collisions are absent from a molecu- 
lar structure and its contactsi 4 i 5 i 6 i 31 Dissipation occurs 
only in the idealized macroscopic reservoirs connected 
by leads to the molecular device. This central role of 
the reservoirs in the process of dissipation is a valid ap- 
proximation, particularly when the applied bias is small 
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and a device operates in a linear regime. When inelastic 
scattering dominates this picture is not applicable. In 
particular, experiments 37 - 38 indicate that inelastic scat- 
tering with lattice vibrations is present at sufficiently 
large bias, causing local heating of contacts and molecu- 
lar devices. Energy transfer to the lattice may also cause 
atomic migration and result in dramatic changes in the 
device characteristics (also they may give phonon replica 
structure in the measured conductance). The modelling 
of a many-electron system out of equilibrium coupled to 
lattice vibrations is a real theoretical challenged 

Electron correlations are also important in molecular 
conductors, for example, Coulomb blockade effects domi- 
nate the transport in quantum dots. Short-range electron 
correlations seems to be relevant in order to get quantita- 
tive description of I/V characteristics in molecular con- 
strictions 42iiii22i In particular it is commonly assumed 
that the energy scales for electron-electron and electron- 
phonon interactions are different and could be treated 
separately. However, the metallic screening of the elec- 
trodes considerably reduces the Coulomb-charging en- 
ergy (from eV to meV scale). In this regime the energy 
scale for the two interactions merge and they need to 
be treated on the same footing posing some additional 
theoretical challenge. 

Other approaches put forward in the literature di- 
rectly look for a homogeneous current-carrying state ei- 
ther based on a a maximum entropy principle^ 3 - or by a 
imposing the current through Lagrange-multipliers^ In 
those approaches it is implicitly assumed that the origin 
of the homogeneous current is independent whether it is 
introduced by reservoirs or by external fields. This is in- 
deed the case for independent electrons but once dissipa- 
tion is built in the system might exhibit a dependence on 
the history of the applied bias, (e. g., possible appearance 
of hysteresis loop in the current versus voltage character- 
istics). 

It is clear that the quality of the TDDFT function- 
als is of crucial importance. In particular, exchange 
and correlation functionals for the non-equilibrium situ- 
ation are required. Time-dependent linear response the- 
ory for dc-steady state has been implemented in Ref. EH 
within TDLDA assuming jcllium-like electrodes (mim- 
icked by complex absorbing/emitting potentials). It has 
been shown that the dc-conductance changes consider- 
ably from the standard Landauer value. Therefore, a 
systematic study of the TDDFT functionals themselves 
is needed. A step beyond standard adiabatic approxima- 
tions and exchange-only potentials is to resort to many- 
body schemes as recently done for the characterization 
of optical properties of semiconductors and insulators. 46 
Another path is to explore in depth the fact that the true 
exchange-correlation potential is current dependent 

An appealing feature of the present approach is that 
electron-electron and electron-ion correlations and dissi- 
pation would in principle be described correctly in two- 
component TDDFTia 

At present we are implementing our propagation 



scheme for real 3D systems 3 ^ within the real-space real- 
time TDDFT code, octopus^ We are also exploring the 
possibility of a semiclassical description of the electron- 
ion coupling in order to avoid the complexity of mul- 
ticomponent DFT and the problems related to mixed 
quantum classical approaches (i.e., Ehrenfest dynamics) 
as they fail to describe the long-term inelastic electron- 
phonon scattering correctly^iSiiS 
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APPENDIX: CONTINUED MATRIX FRACTIONS 



Let us consider the infinite tridiagonal block matrix 



M 



A Bi 
Bi Ai i?2 
B 2 A 2 



(A.l) 



where Ai and Bi are NxN matrices (the argument works 
even for matrices of different matching dimensions). We 
write the inverse matrix AfrT 1 as 



AC 1 = 



Qoo Qoi 
Qio Qi 



(A.2) 



where Qoo is the first NxN block of Mq . It is now 
convenient to introduce the matrix M n obtained by drop- 
ping the first nN lines and nN columns of Mo- Then, in 
terms of the rectangular matrix B n = [_B„,0,0, . . .], we 
have 



1 



1 



" A 


Bi ' 




~A Q 







Mi. 







Mi 



(A.3) 



1 





Mi 



Bi 
Bl 



1 



A B 1 
Bl Mi 



From the above Dyson-like equation it is straightfor- 
ward to obtain Qoo = A^ 1 — A^ 1 BiQio and Qio = 
-Mi 1 BlQ 00 . Solving for Q, 
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1 



100 



A a - BiMi'Bl 



(A.4) 
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One can now proceed along similar lines. We define 



Mf 1 = 



Q21 Qi 



(A.5) 



where Qn is the first N x N block of Mf 1 . Fr om the 
corresponding Dyson equation one finds Qn — A^ 1 — 



A^ 1 B 2 Q 21 and Q 21 = -M^B^Q U . Solving for Q n 

1 



Q 
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(A.6) 



A x - B 2 M 2 ~ 1 B 2 r 
and substituting this result back in Qoo yields 

— L t -. (A.7) 



A0-B1 



Ai-B 2 M~ 1 B. 



Repeating the same steps we end up with the continued 
matrix fraction 



A0-B1- 



A1-B2- 



A2-B3 



-B 2 



-Bi 
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